library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())


#----- Load data 
load(file = 'argData.RData')

# No FEs, htun + tb
m1 <- lm(women_tb_htun ~ mag + female + magAnalysisWoman, data = ar_small)
vcov1 <- vcovHC(m1, 'HC2')
s1  <- sqrt(diag(vcov1))
# No FEs,  tb
m2 <- lm(women_tb ~ mag + female + magAnalysisWoman, data = ar_small)
vcov2 <- vcovHC(m2, 'HC2')
s2  <- sqrt(diag(vcov2))
# No FEs, htun + tb
m3 <- lm(women_htun ~ mag + female + magAnalysisWoman, data = ar_small)
vcov3 <- vcovHC(m3, 'HC2')
s3  <- sqrt(diag(vcov3))

# FEs, htun + tb
data_m1 <- fixest::demean(women_tb_htun + mag + female + magAnalysisWoman
		~ term + party, data = ar_small)
m1a <- lm(women_tb_htun ~ mag + female + magAnalysisWoman - 1, data = data_m1)
vcov1a <- vcovHC(m1a, 'HC2')
s1a  <- sqrt(diag(vcov1a))
# FEs,  tb
data_m2 <- fixest::demean(women_tb + mag + female + magAnalysisWoman
		~ term + party, data = ar_small)
m2a <- lm(women_tb ~ mag + female + magAnalysisWoman - 1, data = data_m2)
vcov2a <- vcovHC(m2a, 'HC2')
s2a  <- sqrt(diag(vcov2a))
# FEs, htun + tb
data_m3 <- fixest::demean(women_htun + mag + female + magAnalysisWoman
		~ term + party, data = ar_small)
m3a <- lm(women_htun ~ mag + female + magAnalysisWoman - 1, data = data_m3)
vcov3a <- vcovHC(m3a, 'HC2')
s3a  <- sqrt(diag(vcov3a))

# Generate Table
stargazer::stargazer(m1, m2, m3, m1a, m2a, m3a, 
	se = list(s1, s2, s3, s1a, s2a, s3a), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes'),
		c('FE by Party', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--1999-2014  Cámara de Diputados de la Nación Argentina",
	label = 't:argentina'
)

### Substantive Effect
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Sims
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims1a <- rmvnorm(n = 1000, mean = coef(m1a), sigma = vcov1a)
sims2a <- rmvnorm(n = 1000, mean = coef(m2a), sigma = vcov2a)
sims3a <- rmvnorm(n = 1000, mean = coef(m3a), sigma = vcov3a)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(cbind(1, Wdata)[2,])) - (sims1 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(cbind(1, Wdata)[2,])) - (sims2 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(cbind(1, Wdata)[2,])) - (sims3 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025))

pred1a_W <- apply((sims1a %*% t(Wdata[2,])) - (sims1a %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2a_W <- apply((sims2a %*% t(Wdata[2,])) - (sims2a %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3a_W <- apply((sims3a %*% t(Wdata[2,])) - (sims3a %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))

# Men
pred1_M <- apply((sims1 %*% t(cbind(1, Mdata)[2,])) - (sims1 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.95, 0.5, 0.05))
pred2_M <- apply((sims2 %*% t(cbind(1, Mdata)[2,])) - (sims2 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.95, 0.5, 0.05))
pred3_M <- apply((sims3 %*% t(cbind(1, Mdata)[2,])) - (sims3 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.95, 0.5, 0.05))

pred1a_M <- apply((sims1a %*% t(Mdata[2,])) - (sims1a %*% t(Mdata[1,])), 2, quantile, c(0.95, 0.5, 0.05))
pred2a_M <- apply((sims2a %*% t(Mdata[2,])) - (sims2a %*% t(Mdata[1,])), 2, quantile, c(0.95, 0.5, 0.05))
pred3a_M <- apply((sims3a %*% t(Mdata[2,])) - (sims3a %*% t(Mdata[1,])), 2, quantile, c(0.95, 0.5, 0.05))

# Building blocks, Table M.40
stargazer::stargazer(c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred1a_W[2,], 2), round(pred2a_W[2,], 2),
	round(pred3a_W[2,], 2)))


stargazer::stargazer(c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred1a_W[3,], 2), ', ', round(pred1a_W[1,], 2), ')'), 
	paste0('(', round(pred2a_W[3,], 2), ', ', round(pred2a_W[1,], 2), ')'), 
	paste0('(', round(pred3a_W[3,], 2), ', ', round(pred3a_W[1,], 2), ')'))
)

stargazer::stargazer(c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred1a_M[2,], 2), round(pred1a_M[2,], 2),
	round(pred3a_M[2,], 2)))

stargazer::stargazer(c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred1a_M[3,], 2), ', ', round(pred1a_M[1,], 2), ')'), 
	paste0('(', round(pred2a_M[3,], 2), ', ', round(pred2a_M[1,], 2), ')'), 
	paste0('(', round(pred3a_M[3,], 2), ', ', round(pred3a_M[1,], 2), ')'))
)

